frailtypack (version 2.2-13)

frailtyPenal for Joint frailty models: Fit Joint Frailty model for recurrent and terminal events using penalized likelihood estimation

Description

Fit a joint frailty model for recurrent and terminal events using a Penalized Likelihood on the hazard function. Left-truncated and right-censored data are allowed. Stratified analysis is not possible. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring. A common frailty term $(\omega_i)$ for the two rates will take into account the heterogeneity in the data, associated with unobserved covariates. The frailty term acts differently for the two rates ( $\omega_i$ for the recurrent rate and $\omega_i^{\alpha}$ for death rate). The covariates could be different for the recurrent rate and death rate.

The joint model for the hazard functions for recurrent event $r_i(.)$ and death $\lambda_i(.)$ is :

$$r_i(t|\omega_i)=\omega_ir_0(t)exp(\bold{\beta_1^{'}Z_i(t)})$$ $$\lambda_i(t|\omega_i)=\omega_i^{\alpha}\lambda_0(t)exp(\bold{\beta_2^{'}Z_i(t)})$$

where $r_0(t)$ (resp. $\lambda_0(t)$) is the recurrent (resp. terminal) event baseline hazard function, $\bold{\beta_1}$ (resp. $\bold{\beta_2}$) the regression coefficient vector, $\bold{Z_i(t)}$ the covariate vector. And where $\omega_i\sim\bold{\Gamma}(\frac{1}{\theta},\frac{1}{\theta})$ is iid.

Arguments

formula

a formula object, with the response on the left of a $\texttildelow$ operator, and the terms on the right. The response must be a survival object as returned by the 'Surv' function like in survival package.

formula.terminalEvent

a formula object, only requires terms on the right to indicate which variables are modelling the terminal event.

data

a 'data.frame' in which to interpret the variables named in the 'formula' and 'formula.terminalEvent'.

Frailty

Logical value. Is model with frailties fitted? If so, variance of frailty parameter is estimated. If not, Cox proportional hazards model is estimated using Penalized Likelihood on the hazard function. The default is FALS

joint

Logical value. Is joint model fitted? If so, 'formula.terminalEvent' is required. The default is FALSE

recurrentAG

Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent covar

cross.validation

Logical value. Is cross validation procedure used for estimating smoothing parameter? If so a search of the smoothing parameter using cross validation is done, with kappa1 as the seed. The cross validation is not implemen

n.knots

integer giving the number of knots to use. Value required. It corresponds to the (n.knots+2) splines functions for the approximation of the hazard or the survival functions. Number of knots must be between 4 and 20.(See Note)

kappa1

positive smoothing parameter. The coefficient kappa of the integral of the squared second derivative of hazard function in the fit (penalized log likelihood). To obtain a good value for kappa1, a solution is to fit the corresponding

kappa2

positive smoothing parameter for the terminal event rate. To obtain a good value for kappa2, a solution is to fit the corresponding Cox model using cross-validation (See cross.validation) with the death indicator as the event of interest. See

maxit

maximum number of iterations for the Marquardt algorithm. Default is 350

Value

• Parameters estimates of a joint frailty model, more generally a 'fraityPenal' object. Methods defined for 'frailtyPenal' objects are provided for print, plot and summary. The following components are included in a 'frailtyPenal' object for Joint frailty models.
• alpha
  the coefficient associated with the frailty parameter in the terminal hazard function
• call
  The code used for fitting the model
• coef
  the coefficients of the linear predictor, which multiply the columns of the model matrix.
• cross.Val
  Logical value. Is cross validation procedure used for estimating the smoothing parameters?
• formula
  the formula part of the code used for the model
• groups
  the maximum number of groups used in the fit
• kappa
  A vector with the smoothing parameters corresponding to each baseline function as components
• lam
  matrix of hazard estimates at x1 times and confidence bands.
• lam2
  the same value as lam for the second stratum
• n
  the number of observations used in the fit.
• n.events
  the number of events observed in the fit
• n.iter
  number of iterations needed to converge
• nvar
  A vector with the number of covariates of each type of hazard function as components
• surv
  matrix of baseline survival estimates at x1 times and confidence bands.
• surv2
  the same value as surv for the the second stratum
• theta
  variance of the frailty parameter $(\bold{Var}(\omega_i))$
• type
  character string specifying the type of censoring. Possible values are "right", "left", "counting", "interval", "interval2". The default is "right" or "counting" depending on whether the 'time2' argument is absent (not interval-censored data) or present (interval-censored data), respectively.
• varH
  the variance matrix of all parameters (theta, the regression coefficients and the spline coefficients).
• varHIH
  the robust estimation of the variance matrix of all parameters (theta, the regression coefficients and the spline coefficients).
• x1
  vector of times where both survival and hazard function for the recurrent events are estimated. By default seq(0,max(time),length=99), where time is the vector of survival times.
• x2
  see x1 value for the terminal events survival and hazard function

synopsis

frailtyPenal(formula, formula.terminalEvent, data, Frailty = FALSE, joint = FALSE, recurrentAG = FALSE, cross.validation = FALSE, n.knots, kappa1, kappa2, maxit = 350)

Details

The estimated parameter are obtained by maximizing the penalized loglikelihood using the robust Marquardt algorithm (Marquardt, 1963) which is a combination between a Newton-Raphson algorithm and a steepest descent algorithm. When frailty parameter is small, numerical problems may arise. To solve this problem, an alternative formula of the penalized log-likelihood is used (see Rondeau, 2003 for further details). Cubic M-splines of order 4 are used for the hazard function, and I-splines (integrated M-splines) are used for the cumulative hazard function. The inverse of the hessian matrix is the variance estimator and to deal with thepositivity constraint of the variance component, a squared transformation is used and the Standard errors are computed by the $\Delta$-method (Knight & Xekalaki, 2000). The smooth parameter can be chosen by maximizing a likelihood cross validation criterion (Joly and other, 1998). The integrations in the full log likelihood werre evaluated using Gaussian quadrature. Laguerre polynomials with 20 points were used to treat the integrations on $[0,\infty[$.

INITIAL VALUES

The splines and the regression coefficients are initialized to 0.5. The program fits an adjusted Cox model to have new initial values for the regression and the splines coefficients. The variance of the frailty term $\theta$ and the coefficient $\alpha$ associated in the death hazard function are initialized to 1. Then, it fits a joint frailty model.

PARAMETERS LIMIT VALUES

As frailtypack is written in Fortran 77 some parameters had to be hard coded in.The default values of these parameters are, with the corresponding variable namein the fortran code between brackets.

maximum number of observations(ndatemax): 30000 maximum number of groups(ngmax): 1000 maximum number of subjects(nsujetmax): 15000 maximum number of parameters (npmax) :50 maximum number of covariates (nvarmax):50 If these parameters are not large enough (an error message will let you know this), you need to reset them in joint.f and recompile.

References

V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P. Soubeyran (2007). Joint frailty models for recurring events and death using maximum penalized likelihood estimation:application on cancer events. Biostatistics, 8,4, 708-721.

V. Rondeau, D Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.

D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.

See Also

summary.jointPenal, print.jointPenal, plot.jointPenal, readmission, terminal, cluster

Examples

### Joint model (recurrent and terminal events) with 2 covariates ###
### on a simulated dataset ###


data(readmission)

## Gap-time ##
modJoint_gap<-frailtyPenal(Surv(time,event)~cluster(id)+sex+as.factor(dukes)
		+as.factor(charlson)+terminal(death),
		formula.terminalEvent=~sex+as.factor(dukes)+as.factor(charlson),
		data=readmission,n.knots=14,kappa1=9550000000,
		kappa2=1410000000000,Frailty=TRUE,joint=TRUE,recurrentAG=FALSE)

## Calendar time ##
modJoint_calendar<-frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)+sex
		+as.factor(dukes)+as.factor(charlson)+terminal(death),
		formula.terminalEvent=~sex+as.factor(dukes)+as.factor(charlson),
		data=readmission,n.knots=10,kappa1=9550000000,
		kappa2=1410000000000,Frailty=TRUE,joint=TRUE,recurrentAG=TRUE)
    
  print(modJoint_gap)
  summary(modJoint_gap)
  plot(modJoint_gap)
  
  print(modJoint_calendar)
  summary(modJoint_calendar)
  plot(modJoint_calendar)

# A model takes around 1 minute to converge #